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Abstract 

The objective of this article is to complete preliminary results from 
[H [13] concerning the time-minimal control of dissipative two-level quan- 
tum systems whose dynamics is governed by Lindblad equations. The 
extremal system is described by a 3D-Hamiltonian depending upon three 
parameters. We combine geometric techniques with numerical simulations 
to deduce the optimal solutions. 
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1 Introduction 

In this article, we consider the time-minimal control analysis of two-level dis- 
sipative quantum systems whose dynamics is governed by Lindblad equation. 
More generally, according to [lOi, the dynamics of a finite-dimensional quantum 
system in contact with a dissipative environment is described by the evolution 
of the density matrix p given by the equation 

i^^=[H„ + Hi,p]+iC{p), (1) 

where Hq is the field-free Hamiltonian of the system. Hi represents the inter- 
action with the control field and C the dissipative part of the equation; [A, B] 
is the commutator of the operators A and B defined by [A, B] — AB — BA. 
Equation ^ is written in units such that h = I. The components of the density 
matrix satisfy the following equations: 

Pnn = -i[HQ + Hi,p]nn - ^ IknPnn + ^ Inkpkk 

k^n k^n 

Pkn = -i[Ho+ Hi,p]kn-TknPkn, k^n (2) 
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where 1 < k < N and 1 < n < N for a iV-level quantum system. The parameters 
7fc„ describe the population relaxation from state k to state n whereas Tkn is 
the dephasing rate of the transition fr'om state k to state n. Note that not 
every positive parameter 7fe„ or Tkn is acceptable from a physical point of view 
since the density matrix p must satisfy particular properties (trace conservation, 
hermitian operator and complete positivity pH^fTH]). 

Particularizing now to the case N = 2, we assume that Hi is of the form 

Hi = —pxEx — fJ-yEy, 

where the operators fix and are proportional to the Pauli matrices ax and 
(Ty in the eigenbasis of Hq. The electric field is the superposition of two linearly 
polarized fields Ex and Ey and we assume that these two fields are in resonance 
with the Bohr frequency E2 — E1. In the RWA approximation, the time evolution 
of p{t) satisfies the following form of the Lindblad equation 
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where E is equal to E 
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and u is the complex Rabi frequency of the laser 



field (the real and imaginary parts of u are the amplitudes of the real fields 
Ex and Ey up to a multiplicative constant). In Eq. ([3]), is the difference of 
energy between the ground and excited states of the system. In the interaction 
representation, Eq. Q becomes 
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The interaction representation means that we have transformed the mixed-state 
p with the unitary transformation U = diag(l, e*"*, e~"^*, 1). Since Tr[p] = 1, 
the density matrix p can be represented by the vector q (x, y, z) where 
X = 29fi[pi2], y — 23[/9i2] and z = P22 ~ Pii and q belongs to the Block ball 
\q\ < 1. Equation (HJ takes the form: 



X = — r.T + U2Z 

y = -ry - uiz 

z = (712 - 721) - (712 + 721)2^ + uiy - U2X 



(5) 



A = (r,7+,7_) is the set of parameters such that 7+ = 712 +721 and 7_ — 
712 ~ 721 and they satisfy the following inequations 2r > 7+ > |7_| derived 
from the Lindblad equation [T^], the Bloch ball [g] < 1 being invariant. The 
control is u = ui + iu2 where ui and U2 are two real functions. We can write 
the control field u = |u|e"^ where |u| < M and up to a rescaling of the time and 
dissipative parameters we can assume that |u| < 1. 

We consider the time- minimal transfer problem from a state go to a state qi . 
Hence, we have to analyze a time-minimal control problem for a bilinear system 
of the form: 



q^Foiq) 



2 

E 



u^Fi{q), \u\ < 1, 
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where the drift term _Fo depends upon three parameters. This problem is a 
very difficult problem whose analysis requires advanced mathematical tools from 
geometric control theory and numerical simulations. 

Such analysis is motivated by physical reasons. It is a fundamental model in 
quantum control and more general problems can be handled by coupling such 
systems. Numerous optimal control results exist in the conservative case e.g. 
[5], but only partial ones for this problem: a pioneering work [T3] assuming u 
real and a second one [1] for u complex but restricted to 7_ = 0. 

The objective of this article is double. First of all, we introduce all the proper 
geometric tools to analyze the problem using Pontryagin maximum principle. 
Secondly, our aim is to make a complete study for every generic parameter 
in A, combining mathematical reasonings and numerical simulations based on 
shooting techniques and including computations of conjugate points to test op- 
timality. Based on the Cotcot code ^ , they can be used in practice to compute 
the true optimal control, once the physical parameters are identified. 

The organization of this article is the following. In section O we complete 
the classification of the time-minimal synthesis of |13| corresponding to the case 
where u is real but introducing more general tools to handle the problem. It 
corresponds to a time-minimal control problem of a two-dimensional bilinear 
system in the single input case. The optimal synthesis for a fixed initial point 
can be constructed by gluing together local optimal syntheses. We can also 
make estimates of switching points by lifting the system on a semi-direct Lie 
group. This classification is physically relevant to analyze the 3D-case because, 
using the symmetry of revolution of the problem, it gives the time-minimal syn- 
thesis for initial points qo (0,0, ±1). This geometric property is explained 
in section [3l Moreover, using spherical coordinates the system can be viewed 
as a system on a two-sphere of revolution coupled with the evolution of the 
distance to the origin, which represents the purity of the system. According 
to the maximum principle, smooth extremals are solutions of the Hamiltonian 
vector field ll where H ^ Hq + {Hi + Hi = {p,F,{q)), and the con- 

trol components are given by Ui = Hi/{Hf -\- i = 1,2. Non smooth 

extremals can be constructed by connecting smooth subarcs of the switching 
surface E: Hi = H2 = 0. A contribution of this article in section [3] is to classify 
the possible connections. We proved that every non smooth extremal is either a 
solution of the 2D-single input system, assuming u real, or occurs when meeting 
the equatorial plane of the Bloch ball. In the second case, the switching can 
be handled numerically using an integrator with an adaptative step. In the 
same section, we combine analytical and numerical analysis to determine the 
extremals and compute conjugate points. This completes the analysis from [3] 
in the integrable case. The physical interpretation is presented as a conclusion. 

2 The 2D-case 

Following [T^ , a first step in the analysis is to consider the following reduced sys- 
tem. Assuming u real, the x-coordinate is not controllable and we can consider 
the planar single-input system: 

y = -Fy-Miz 

z ^ 7_ - 7+z -I- uiy, < 1. 
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It gives the time-optimal analysis of the control problem when the initial state 
q{0) = (2/(0), z(0)) is a pure state on the z-axis, that is q{0) = (0,±1). We 
proceed as follows to make the analysis. 



2.1 Symmetry group 

A discrete symmetry group is associated to reflections with respect to the differ- 
ent axes. More precisely, if w ^ —z then one gets the same system changing ui 
into —ui and 7_ into — 7~. If w = —y then one gets the same system changing 
ui into —Til. In particular, concerning the time-minimal control problem, this 
amounts to exchange the trajectories (T+ and ct_ corresponding respectively to 
ui = 1 and ui — —1. Also, according to this property, the time- minimal synthe- 
sis is symmetric with respect to the z-axis. This is connected to the symmetry 
of revolution of the whole system which is explained later. 



2.2 The feedback classification 

A preliminary step in our analysis is to consider the feedback classification 
problem |3J. The system is written in a more compact form as follows: 

q = Fiq)+uG{q) 

where F and G are affine vector fields. To make the feedback classification, we 
relax the control bound \u\ < 1. The geometric invariants are related to the 
sets: 

• S = {q, det(G, [F, G]) = 0} where are located the singular trajectories. 

• G — {q, det(F, G) — 0} corresponding to the set of points where F and G 
are coUinear. 

They are obtained by straightforward computations of Lie brackets: 
Hence, S is given by: 

j;[2(r-7+)z + 7_] =0 

and if 7+ 7^ F, the singular set is defined by the two lines y = and z = 
7_/[2(7+ — r)]. The collinear set is defined by: 

j+z^+Ty^^j^z = 

and is a closed curve formed by the union of two arcs of parabola, containing 
(0,0) and (2:1,0), with zi = 7-/7+, which is the equihbrium state of the free 
motion. More generally, C contains the equilibrium points for the dynamics 
with constant control uq, since F + uqG = 0. In particular, the equilibrium 
point when mi = 1 is given by 

Ci : y = r-^— ^ = -ry. 

1 + 7+r 
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Figure 1: Diagram of the sets S and C in solid lines for 7_ = —0.2, 7+ = 0.4 
and r = 1. The equation of the horizontal dashed line is z = 7_/2(7_|_ — F). 



The collinear set C shrinks into a point when 7_ — 0. Computing the intersec- 
tion of C with the singular line z — 7_/[2(7+ — F)], one gets: 



1- 



4(7+ - F)2 



(7+-2F). 



If 7_ ^ 0, since 2F > 7_|_, one deduces that the intersection is empty except in 
the case where 7+ = 2F. An important consequence is to simplify the classifi- 
cation of the optimal syntheses. We represent on Fig. [T]the sets S and C for a 
situation with 7_ < and 7+ — F < 0. Another feedback invariant is the op- 
timality status of singular trajectories. Using the generalized Legendre-Clebsch 
condition, it is splitted into fast and slow directions. In the 2D-case, it is tested 
by Lie brackets configurations and can be computed by introducing: 

D = det(G, [[G, F],G]), D" = det(G, F). 

The trajectory is time-optimal in the so-called hyperbolic case DD" > and 
time-maximal in the so-called elliptic case DD" < 0. Computing Lie brackets 
of length 3, we have: 



and 



[[G,F],F] 
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For the singular direction y = 0, we get: 

DD" = 2z\j+-rh+{z 



2(7+- 
7-2(74 



-F)z 
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Near the origin, the sign is always positive if 7_ 7^ 0. If 7_ = 0, the sign is 
given by (7+ — F). For the singular direction z = 7_/[2(7-|_ — F)], we have: 

DD" = - 2r) - iTy'ii+ - rn 
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Hence, near the origin, the optimahty is given by the sign of (7+ — 2r)(7+ — F). 
Moreover, since T > 7+/2, one gets that DD" > if 7+ - T < and DD" < 

if 7+ - r > 0. 

2.3 Computation of a normal form 

In order to analyze the time-optimal control, we compute a global normal form 
up to a polar singularity for the action of the feedback group. A first step 
is to linearize the vector field G since it is connected to the evaluation of the 
switching times. One further normalization is to straight the horizontal singular 
line: z = 7_/[2(7+ — F)]. Using polar coordinates: 

y = p cos (j), z = p sin cj), 

one gets: 

p = 7_ sin (?!) + p[-r + (r - 7+) sin^ (/)] 

• 7- cos (/) sm(2^ 

^ = — - — + (r-7+) — - — +ui. 
P ^ 

If we use the coordinates x = p^/2 and z, the system becomes: 

X = -2ra; + 7_^; + 2;^(r-7+) 
z = 7_ — 7+^; + u\\j2x — z^. 

Making a feedback transformation of the form U\ — » (iu\ where /3 is a function 
of X and z, we can consider the system: 

X = -2ra; + 7_^ + 2;^(r-7+) 
z = 7_-7+^; + «i. 

If we set ^; = Z + zq where zq = 7_/[2(7+ — F)], we obtain the system: 

7_(7+-2r) 

In this simplified model where the control is resealed by the positive parameter 
P, we keep most of the information about the initial system. In particular, all 
the feedback invariants are preserved: the collinear set corresponds to a; = 
and the singular set is identified to z = with its optimality status, while we 
have wiped out the singular line y = 0. Due to the feedback transformation, we 
lose however the singularities of the vector fields F + G and F — G and also the 
saturation phenomenon of the singular control. 

For the simplified model, the adjoint system takes the form: 

= 2Tpx 

= -2zp:j;(F -7+) +p^7_, 
and can be easily integrated to compute the time-minimal synthesis. 
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2.4 The saturation phenomenon 

One interesting property which is not captured by the normal form is when the 
singular control is saturating along the horizontal singular line z = 7_/[2(7+ — 
r)]. Introducing D' = det(G', [[G, F], i^]), we get on the singular hne: D' = 
2/7- (2r — 7+). The singular control is given by: 

D^_ 7-(7+-2r) 
D- 2y(r-7+) ' 

and saturation occurs when \us\ — 1. Observe that if 7-(2r — 7+) ^ then the 
singular control is never admissible when ?/ = 0. 



2.5 The switching function 

For the true system, the switching function is more intricate but can be still 
analyzed using geometric and numerical techniques. The switching function 
$ is given by = p{t)G{q{t)) and switching occurs when <I>(i) = 0. By 
construction, G is tangent to the circle , hence is rotating when we follow an 
arc curve (T+ or (t_ . The dynamics of p is given by the adjoint equation: 

.dF dG. 
P = -P(^ + ^^)^ " = ±1- 

The corresponding dynamics is linear and p can be either oscillating if the 
eigenvalues are complex or non oscillating if they are real. 

An equivalent but more geometric test is the use of the standard (^-function 
introduced in [7] and defined as follows. Let v be the tangent vector solution of 
the variational equation: 

• dG. 

whose dynamics is similar to the one of the adjoint vector. Let ii < ^2 be two 
consecutive switching times on an arc cr+ or cr_. One can set ti = and t2 = t. 
By definition, we have: 

p{0)G{qm^p{t)G{q{t))^0. 

We denote by «(•) the solution of the variational equation such that v{t) — 
G{q{t)) and where this equation is integrated backwards from time t to time 0. 
By construction p{0)v{0) = and we deduce that at time 0, p{Q) is orthogonal to 
G{q{Q)) and to w(0). Therefore, v{0) and G(g(0)) are collinear; 9{t) is defined as 
the angle between G{q{Q)) and v{Q) measured counterclockwise. One deduces 
that switching occurs when 9{t) = [tt]. In the analytic case, 9{t) can be 
computed using Lie brackets. Indeed, for m = e, e = ±1, we have by definition 

v{0) = e-'^'^^^+'^^G{q{t)), 

and in the analytic case, the ad-formulae gives : 

^(0) = E ^ad"(i^ + eG)G{q{t)). 



n>0 
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Here, to make the computation explicit, we take advantage of the fact that we 

can lift our bilinear system into an invariant system onto the semi-direct Lie 
group GL(2,R) xg identified to the set of matrices of GL{3,'M.): 



1 
9 V 



geGL{2,M.), v€ 
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acting on the subspace of vectors in M"* ^ 

Lie brackets computations are defined as follows. We set: 

F{q)=Aq + a, G{q) = Bq, 

and F,G are identified to {A, a), {B,0) in the Lie algebra gl{2,R) x M?. The 
Lie brackets computations on the semi-direct product are defined by: 

[{A',a'),{B',b')] = {[A',BlA'b'-B'a'). 

We now compute exp[— tad(F + eG)]. The first step consists in determining 
exp[—tad{A + eB)] which amounts to compute a,d{A + sB). We write gl{2,R) = 
c © sl{2, R) where c is the center 

1 
1 

We choose the following basis of s{(2,R): 

The matrix A is decomposed into: 

. f -T 0\ f \ \ / s 

and hence A = -(r + 7+)/2 and s = (7+ - T)/2. In the basis {B, C, D), 
ad(j4 -|- eB) is represented by the matrix: 

-2s 
-2s 2£ 
-2e 

The characteristic polynomial is -P(A) = — A(A^ + 4(e^ — s^)) and the eigenval- 
ues are A = and Xi = ±2\/s^ — e^, i = 1,2; Ai and A2 are distinct and real if 
I7+ — r| > 2 and we note Ai = 2vfc^-s^, A2 = — Ai; Ai and A2 are distinct and 
imaginary if |7_|_ — Fj < 2 and we note Ai = 2i\/£^ — s^, A2 = — Ai. To compute 
g-fad(A+£B)^ we must distinct two cases. 

Real case: In the basis B, C, D, the eigenvectors corresponding to {0, Ai, A2} 
are respectively: Vq =* (e, 0,s), vi =* (2s, — Ai,2e) and V2 =* (2s,— A2,2e). 
Therefore, in this eigenvector basis, exp[—ta,d{A + eB)] is the diagonal matrix: 
diag(l, e~^^*, e~^'^*). To compute exp[— iad(^ + sB)]B, we use the decomposi- 
tion 

B = avo + Pvi + JV2, 
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with: 

" ~ e2 _ ^2 ' - 2(A2 - Ai)(e2 - s^) ' T - 2(Ai - A2)(e2 - s^) ' 
Hence one gets: 

g-tad(A+si3)^ = ai;o + /3e"^^*wi + 7e-^^*V2. 

To test the coUinearity at go, we compute 

det(B(go),e-*^'i(^+^^)i?(9o)) = 

where the determinant is equal to 

(z2 - + 2£(/3e-^i* + 7e-^=*)) + 2yo^o(Ai/3e-^^* + AaTe'^^*). 

Using the fact that this last expression has at most two zeros, one being for 
i = 0, it is straightforward to check that for qq =* (0, 1), this expression only 
vanishes at t = 0. This proves the result numerically checked in [13 that there is 
at most one switching. This analysis can be generalized to any initial condition. 
Imaginary case: In this case, we note Ai = i9 the eigenvalue associated to the 
eigenvector *(2s, — i6', 2e). We consider the real part vi =* (2s,0,2e) and the 
imaginary part V2 =* (0,-0,0). In the basis vq (e, 0,s), ui, W2, ad(A + ei?) 
takes the normal form: 



diag(0. 

Hence, we have in this basis: 



9 
-9 



cos(e't) - sin(6't) 
sin(6't) cos(e'i) 



We decompose B in the same basis: 

B — avo + [3vi + 7^2, 

where 

Hence, we get: 

^_tad(A+eB)^ = avQ + l3[cos{9t)vi + sin(e'i)w2]. 
Computing we obtain: 
det(Bgo,e^*''^(^+^-^)B(qo)) = (^o " vDi^s + 2£/3cos(6't)) + 2/361 sin(e't)yo^o 

which vanishes for two values of t in [0, 27r/0[ if j/o 7^ ^o- These two values 
coincide for =* (0,1). Hence the switchings occur periodically with a period 
of 211/9 which confirms the numerical simulations of [IS]. 
Case 7_ ^ 0: The Lie bracket is given by: 

[{A',a'),{B',b')]^{[A',B'],A'b' ~B'a'). 
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We note (ei, 62) the M^-canonical basis and a = 7-62. Wc have: 

a.d{A + eB,a)- {B,0) = {[A + eB,B],-Ba) 
= i[A,B],-Ba), 

a<f{A + eB,a)-{B,0) = [{A + eB,a),{[A, B], -Ba)] 

^ {ad^iA + eB, B), -{A + eB)Ba - [A, B]a). 

More generaUy, one gets: 

a.A^{A + eB, a) ■ (B, 0) = {ad''{A + eB) ■ B, Vk) 

where Vk is given by the recurrence relation: 

Vk = -ad''~^(A + eB) ■ Ba+{A + eB)vk-i. 

The computation is intricate but simphfies if 7+ — T since in this case 
ad'' {A + eB) ■ B = for A; > 1. Numerical simulations have to be used to 
compute the switchings sequence. 

Generalization: This technique can be generalized to the time-minimal con- 
trol problem in the full control case, replacing the control domain |m| < 1 by 

\UI\,\U2\ < 1. 

2.6 The time- minimal synthesis 

We use [3] as general reference on time-minimal synthesis, see also [7]- The initial 
condition is fixed to qq (0, 1) and we consider the problem of constructing 
the time-minimal synthesis from this initial point. This amounts to compute 
two objects: 

• The switching locus ^(qo) of optimal trajectories which is deduced from 
the switching locus of extremal trajectories. 

• The cut locus C{qo) which is formed by the set of points where a minimizer 
ceases to be optimal. 

In order to achieve this task, we must glue together local time minimal syntheses 
which are classified. To be more precise, we recover the case (d) of ^13j, the 
gluing being indicated on Fig. [5]on which we have represented the local extremal 
classifications of [5] which are crucial to deduce the optimal syntheses. In this 
case, the cut locus is a segment of the z— axis and its birth is located at the 
initial point (0, 1) which is a consequence of the elliptic situation. The switching 
locus is the union of a fast singular trajectory, corresponding to an hyperbolic 
point and a curve Si (go) corresponding to parabolic points (for the terminology 
see [3]). Note also the importance of the tangential point where arcs (t+ and 
(T_ are tangent leading to the fish-shaped accessibility set ^+(go) represented 
on Fig. [21 This set is not closed since the arc starting from (0, zi) is not 
in A'^{qo). We next list the micro- local situations we need to construct the 
synthesis. 

• Ordinary switching points: The local synthesis is given by a-a+ or 
(T_i_CT_. The two cases arc distinguished using for instance the clock form 
u = pdq with (p, G) — and (p, F) — I which is also useful to get more 
global results. 
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Figure 2: Cut and switching loci for the case 7_ < 0. The sets C and S are 
represented in dashed Hnes, the switching locus in dotted lines. 
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Figure 3: Poisson shape of the accessibihty set. 
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Figure 4: Structure of the extremals for the hypcrbohc case. 




Cut locus 



Figure 5: EUiptic case. 

• Hyperbolic case: Existence of a fast singular admissible trajectory. The 
optimal synthesis is of the form a±as(y± where CTs is a singular arc (see 
Fig.Sl). 

• Elliptic case: Existence of a slow admissible singular trajectory. An 
optimal arc is bang-bang with at most one switching. Not every extremal 
trajectory is optimal and we have birth of a cut locus (see Fig. 

• Parabolic point: It corresponds to a non-admissible singular direction. 
Every extremal curve is bang-bang with at most two switchings. In our 
case, the initial point is fixed and the switching locus starts with the 
intersection of (t_ with the singular line (see Fig. [6]). 

• Saturating case: 

A fast singular trajectory is saturating at a point M : birth of a switching 
curve at M (see Fig. [7]). 

• A C n S" ^ case: 

A fast singular trajectory meets the set C and becomes slow. 




Singular direction 



Figure 6: Parabolic case. 
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Switching locus 



Figure 7: Saturating case 





/ Collinear set 




■m — Singular line 



Figure 8: A C n S" 7^ case 



2.7 Classification of the optimal syntheses 

We describe the different time-optimal syntheses in the single-input case. With- 
out loss of generality, we restrict the study to the initial point qo =* (0, 1). The 
classification is done with respect to the relative positions of the feedback in- 
variants C and S and to the optimal status of singular extremals which are fast 
or slow according to the values of F, 7+ and 7_. 

For 7_ = 0, the set C is restricted to the origin and we have two cases 
according to the sign of F — 7+ . Note that the form of the extremals (T+ and 
CT_ starting from go depends on the sign of |F — 7+I — 2. Two cases for F > 7+ 
are presented in [T^]. We complete this study with the optimal synthesis for 
F < 7+ and |F - 7+I < 2 displayed in Fig. Ek- 

For 7_ 7^ 0, we distinguish four cases according to the signs of 7_ and F — 7_|_. 
One case (F > 7+ and 7_ < 0) is treated in [13]. We consider here three types 
of optimal synthesis represented in Fig. O), O; and [HJi. Note that in a same 
class of synthesis the reachable set from the initial point qq depends on the 
dissipative parameters which can modify the structure of the synthesis. The 
last case 7_ > and 7+ > F can be deduced from the case 7_ > and 7+ < F 
since the horizontal singular line plays no role in both cases. The synthesis of 
Fig. is very similar to the one of Fig. [5] except the fact that a part of the 
horizontal singular line is admissible. The switching locus has been computed 
numerically using the switching function <i>. 

The role of the parameter 7„ is clearly illustrated in Figs. UK and [Hfc. The 
case 7_ = is a degenerate case where the set C shrinks into a point. The vari- 
ation of 7_ induces a bifurcation of the control system leading to new structures 
of the optimal synthesis. For 7_ 7^ 0, the set C is the union of two branches 
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Figure 9: Optimal syntheses for (a) (F = 1.1, 7+ = 1.6, 7_ = 0), (b) {T = 4, 
7+ = 1.5, 7_ = 0.5), (c) (r = 4, 7+ = 6.5, 7_ = -1.5) and (d) (T ^ 1, 
7+ — 0.5, 7_ = —0.1). Solid and dashed vertical and horizontal lines correspond 
respectively to fast and slow singular lines. The set C is represented in dashed 
lines. The switching locus is plotted in dotted line. In (d), only the admissible 
singular horizontal line is represented in solid line. 

of parabola. The optimal status of the vertical singular line changes when this 
line crosses the set C in Fig. [SJ;. 

3 The bi-input case 
3.1 Geometric analysis 

The system is written in short in Cartesian coordinates as follows: 

q = Fo{q) + UiFi{q) + U2F2{q), \u\ < 1. 

Introducing the Hamiltonians Hi — {p, Fi), z = 0, 1, 2, the pseudo-Hamiltonian 
associated to the time-optimal control problem is: 

2 

1=1 

where pa < 0. The time-optimal control is given outside the switching surface 
E: Hi = H2 = 0, by Ui — Hi/ \J H\ -I- H\ , i = 1,2, with the corresponding true 
Hamiltonian: 

Hr^Ho + ^Hl+Hl, 

whose solutions (outside S) are smooth and are called extremals of order 0. 
More general non smooth extremals can be obtained by connecting such arcs 
through S. 
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To make the geometric analysis and to highlight the symmetry of revolution, 
the system is written using the spherical coordinates: 



X — p sm (j) cos 0, y — p sm (p sm &, z — p cos qy 
and a feedback transformation: 

vi — ui cos 6 + U2 sin 9, V2 — — ui sin 6 + U2 cos 9. 
We obtain the system: 

p = 7_ cos (j) — p(7+ cos^ + r sin^ 0) (6a) 

■ 7_sin(/) sin(20) 

= + — (7+-r)+W2 (6b) 

p 2 

9 = — cot(/)i;i. (6c) 
Hence, one deduces that the true Hamiltonian is: 

Hr = [7- cos0-p(7+ cos2 0+rsin2 0)]Pp+P0[- + ^'''^^^'^^ (7+ -r)] + + ( 
From this, we deduce the following lemma: 

Lemma 1. (i)- The angle 9 is a cyclic variable and pg is a first integral (sym- 
metry of revolution). 

(a)- For 7_ — 0, using the coordinate r ~ In p. the Hamiltonian takes the form: 



Hr = -(7+ cos^ (/) + r sin^ <l))pr + sin(20)(7+ - T)p^ + ^pI+pI cot^ 0. (7) 

Hence, r is an additional cyclic variable and pr is a first integral. The system 
is thus Liouville integrable. 

As a consequence, we can deduce two properties. First of all, the z-axis is 
an axis of revolution and the state go =* (0,0, 1) is a pole. This means that by 
making a rotation around {Oz) of the extremal synthesis for the 2D-system, we 
generate the extremal synthesis for the 3D-system. 

More generally, we have for 7_ = a system on the two-sphere of revolution 
described by Eqs. (I6bp and ((6c|) coupled with the one dimensional system ((6a| 
describing the evolution of the physical variable p corresponding to the purity of 
the system. Moreover, the system is invariant for the transformation (/) 1— > tt — 
which is associated to a reflexional symmetry with respect to the equator for the 
system ^ restricted to the two-sphere of revolution. This property is crucial 
in the analysis of the integrable case. 

If 7_ ^ then the situation is more intricate. The extremals solutions of 
order satisfy the equations which arc singular for p — Q: 

p = 7_ cos (f> — p(7+ cos^ 4> + T sm^ 0) 

4> - _I^ + 5i^lM(,^_r) + f^ (8) 

p 2 Q 

■ pe cot^ (j) 
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and 



Vp = (7+ cos ^ + r sm (pjpp — — 

P4> = b- sin 4> + p{^-l+) sin(2(^)]pp - [- - cos ch- + (7+ - T) cos(20)]p0 + ^^'^"a'^ 
P0 = 0, 



where Q = ^/pfcot^ 



3.2 Regularity analysis 

The smooth extremal curves solutions of Hr are not the only extremals because 
more complicated behaviors are due to the existence of the switching surface E: 
H\ = H2 = 0. Hence, in order to get singularity results, wc must analyze the 
possible connections of two smooth extremals crossing E to generate a piecewise 
smooth extremal. This can also generate complex singularities of the Fuller 
type, where the switching times accumulate. In our problem, the situation is 
less complex because of the symmetry of revolution. The aim of this section is 
to make the singularity analysis of the extremals near E. 

The structure of optimal trajectories is described by the following proposi- 
tion. 

Proposition 1. Every optimal trajectory is: 

• Either an extremal trajectory with pe = contained in a meridian plane 
and time-optimal solution of the 2D-system, where u = ('Ui,0). 

• Or subarcs solutions of Hr, where pg ^ with possible connections in the 
equator plane for which (j) = -k/2. 

Proof. The first assertion is clear. If = then extremals are such that ^ = 

and up to a rotation around the z-axis, they correspond to solutions of the 2D- 
system. The switching surface E is defined by: pe cot (j) = p^ = 0. We cannot 
connect an extremal with pe ^ to an extremal where pe = since at the 
connection the adjoint vector has to be continuous. Hence, the only remaining 
possibility is to connect subarcs of Hr with p$ ^ a,t a, point of E leading to 
the conditions Ptj, = and 4> = 1^/2. □ 

Further work is necessary to analyze the behaviors of such extremals near 

E. 

Normal form: A first step in the analysis is to construct a normal form. Taking 
the system in spherical coordinates and setting tp = 7r/2 — (j), the approximation 
is: 

p = 7-V'-p[r + (7+-r)v'] 

^ = ^(l-V;V2)-V(7+-r)-t;2 

e = -ipvi, 
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with the corresponding Hamiltonian: 



Hr =Pp[7-^-p(r+(7+-r)V'2)]+p^[^(i-vV2)-V'(7+-r)]+y^p| + 

" (9) 

Proposition 2. Near 'ip = 0, = 0, we have two distinct cases for optimal 
trajectories: 

• If J- = 0, for the 2D-system, the line tp = is a singular trajectory with 
admissible zero control if j+ — T ^ 0. It is slow if (7+ — F) > and fast 
if (7-1- — r) < 0. Hence, for this system, we get only extremal trajectories 
through E in the case (7+ — F) < 0, where tp is of order t and p^ of order 
t^ . They are the only non-smooth optimal trajectories passing through S. 

• If')-^ 0, for the 2D-system, the settp = p^ = becomes a set of ordinary 

switching points where -0 and p^ are of order t. Moreover, connections 

for extremals of are eventually possible, depending upon the set of 
parameters and initial conditions. 

Proof. For the normal form, the adjoint system is: 

p, = p,(r + (7+-r)v2) + 57-(i-^) 

'Y—tp 

Pi' = -Ppil- - ^■^p{i+ - r)) +Pi,{—^ + (7+ - r)) + vipe- 

(10) 

In order to make the evaluation of smooth arcs reaching or departing from E, 
the technique is simple: a solution of the form 'ip(t) — at + o{t), p^{i) — bt + o{t) 
is plug in the equations to determine the coefficients. From the equations, we 
observe that the contacts with S differ in the case 7_ = from the case 7_ ^ 
that we discuss separately. 

First of all, we consider the case 7_ = 0: pe — 0, tp = is an admissible 
singular direction (with zero control) which can be slow if (7+ — F) > or fast 
if (7+ — F) < 0. In the first case, there is no admissible extremal through E 
while it is possible if 7+ — F < 0. If we compute the different orders, we have 
that tp is of order t, p^ is of order t^ while Pp has to be non zero if pe = 0. If we 
consider extremals with p^ ^ 0, we can conclude with the orders alone. Indeed 
the Hamiltonian \s = £, £ = Q,! and in both cases, we have: 

-Ppp{l+ - F)V;2 - p^V(7+ - r) + ^Jpl+Pir = 0. 

The conclusion using orders is then straightforward. For instance, if ip and p^ 
are of order one, this gives p^ = pg-tp = which is impossible. The other cases 
are similar. 

In the case 7_ 7^ 0, the analysis is more intricate and we must analyze the 
equations. We introduce the Hamiltonians: 

Hi = -p0ip, H2=p^. 
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Differentiating Hi and H2 with respect to t, one gets: 

Hi = {Hi,Ho} + V2{Hi,H2} 
H2 = {H2,Ho} + vi{H2,Hi} 

and at a point of S, we obtain the relations: 

Hi = -pe{l- - V2), H2 = 1-Pp - viPe- 

In order to analyze the singularity, we use a polar blowing up: 

Hi — r cos a, H2 — r sin a, 

and we get: 

. pe cos a . 

r = 7_ h Pp sm a\ 

P 

1 pe7-sinQ; 

a = -7_PpC0SQ;H pg\. 

r p 

Hence, the extremals crossing S are given by solving a — while the sign of r 
is given by the first equation above. 

Depending upon the parameters and the initial conditions on {pp,p), the 
equation a — can have at most two distinct solutions on (0, 27r), while in the 
case pg = 0, we get an ordinary switching point for the single-input system. 
The assertion [2] is proved. □ 



3.3 Geometric analysis and numerical solution 

We first analyze the integrable case 7_ = 0. We only present a summary of the 
result of [4j in order to be generalized to the case 7_ 7^ 0. 



3.3.1 The case 7 

The system ([7]) is associated to a system on the two-sphere of revolution of the 
form: 

2 

q = Go{q) + ^UiGi{q). 

i=l 

It defines a Zermelo navigation problem on the two-sphere of revolution where 
the drift term Go represents the current: 

sin(20) d 
Go-^^(7+-r)-, (11) 

and Gi — G2 — — cot (j)-^ form a frame for the metric g = dcfP + tan^ ipdd^ 
which is singular at the equator = 7r/2. The drift can be compensated by a 
feedback with |u| < 1 if I7+ — r| < 2. This leads to the following discussion. 

Case I7+ — r| < 2: In this case, the system reduced to the two-sphere defines 
a Finsler geometric for which the extremals are a deformation of the extremals 
of g = d(fP' + tan^ (j)d9'^ . The main problem properties are described in the next 
proposition. 
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Figure 10: Extremal trajectories for T = 2.5, 7+ = 2 and 7_ = 0. Other 
parameters are taken to be p^{0) = —1 and 2.33, 0(0) = 7r/4, pp = 1 and 
pe = 2. Dashed Unes represent the equator and the antipodal parallel located 
at = 37r/4. 




Figure 11: Extremal trajectories for F — 4.5, 7+ = 2 and 7_ = 0. Dashed 
lines represent the equator and the locus of the fixed points of the dynamics. 
The solid line corresponds to the antipodal parallel. Numerical values of the 
parameters are taken to be (j){0) = 2tt/5, pe = 8 and Pp{0) = 0.25. The different 
initial values of p^ are -50, -10, 0, 2.637, 3, 5, 10 and 50. 



Proposition 3. If for fixed {pr,p0), the level set of Hr = e (e — 0,1) is 
compact without singular point and has a central symmetry with respect to {(f> = 
7r/2,p0 = 0) then it contains a periodic trajectory ((j),p,j,) of period T and if 
p^(0) are distinct, we have two distinct extremal curves q~^{t), q~ (t) starting 
from the same point and intersecting with the same length T/2 at a point such 
that (f){T/2) = TT - 0(0) (see Fig. [70j) . 

Case I7+ — F| > 2: We have two types of extremals characterized by their 
projection on the two-sphere: those occurring in a band near the equator and 
described by proposition [3] and those crossing a band near (p = n/i and with 
asymptotic properties of proposition |H 

Proposition 4. // |F — 7+I > 2 then we have extremal trajectories such that 
0—^0, \ptj,\ +00 when t —>■ +00 while 9—^0. 

Both behaviors are represented on Fig. [TTl 
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3.3.2 The case 7 7^ 

We present numerical results about the behavior of extremal solutions of order 
and conjugate point analysis. 

Extremal trajectories: 

We begin by analyzing the structure of extremal trajectories. The description 
is based on a direct integration of the system ([H]). We observe two different 
asymptotic behaviors corresponding to stationary points of the dynamics which 
are described by the following results. 

Proposition 5. In the case denoted (a) where \p,f,{t)\ — > +00 when t +00, 
the asymptotic stationary points {pf,(f>f,6f) of the dynamics are given by pf — 
|7_|Vl + r7(l+7+r) and4>f = arctan(l/r) i/7_ > or0/ = 7r-arctan(l/r) 
ifl- < 0. 

Proof. We assume that +00 as t — > +00 and that cot((/)) remains 

finite in this limit. One deduces from the system ^ that {pf,(f>f) satisfy the 
following equations: 

7_ cos(/)/ = Pf{"i+ cos^ 0/ + Tsiv? cf)f) 
7_ sin 4> f 



Pf 



= (7+ — r) cos 0/ sin + e, 



where e = ±1 according to the sign of p^. The quotient of the two equations 
leads to 

(7+ — r) COS0/ sin^/ + £ = tan(/)/(7+ cos^ 4>S + Tsin^ (/)/) 
which simplifies into 

tan0/ = -. 

Using the fact that 0/ G]0, 7r[ and 7_ cos ipf > 0, one arrives to (j)f — arctan(l/r) 
if 7_ > and <j)f = tt — arctan(l/r) if 7_ < 0. From the equation 

7_ cos 4>f = pf{'-f+ cos^ 0/ + r sin^ 0/ ) , 

one finally obtains that 

7_^/TTT2 

□ 

Proposition 6. In the case denoted (h) where limt^+00 </>(i) = or tt, the 

asymptotic limit of the dynamics is characterized by pf = |7_|/7+ and 0/ = 
i/ 7_ > or (pf — t: if < 0. 

Proof. Using the relation 

7_ cos 0/ = Pfi'^+ cos^ 0/ + r sin^ 0/ ) , 
one deduces that 7_ cos0/ > and that pf = I7-I/7+ if 0/ = or tt. □ 
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We have numerically checked that if |r — 7+I > 2 then only the case (a) 
is encountered whereas if |r — 7+I < 2, the extremals are described by cases 
(a) and (b) . One particularity of the case (a) is the fact that the limit of the 
dynamics only depends on T and on the sign of 7_ and not on (f>{0) or 7+. The 
structure of the extremals is also simple in case (b) since the limit of is or 
TT independently of the values of F, 7+ or 7_. The different behaviors of the 
extremals are illustrated in Fig. [T^lfor the case |F — 7+I > 2 and in Fig. [HI for 
the case |r — 7+I < 2. The corresponding optimal control fields vi and V2 are 
represented in Fig. [T3]for the case (a) and in Fig [TS] for the case (b). In Fig. 
[T3l note that the control vi tends to whereas V2 is close to -1 for t sufficiently 
large. This is due to the fact that \p^\ — > +00 when t — > +00 and can be easily 
checked from the definition of vi and V2- We observe a similar behavior for the 
case (b) in Fig. [TS] The control field vi acquires here a bang-bang structure 
which is related to the unbounded and oscillatory behavior of ^^(i) (see Fig. [TB|) . 

Conjugate points: 

The Cotcot code is used to evaluate the conjugate points. This occurs only 
in case (b) and the numerical simulations give that the first conjugate points 
appear before an uniform number of oscillations of the (j) variable. This phe- 
nomenon is represented on Fig. 1161 Cutting the trajectory at the first conjugate 
point avoids such a behavior. Note that due to the symmetry of revolution, the 
global optimality is lost for 6 < tt. 




Figure 12: Extremal trajectories for F = 4.5, 7+ = 2 and 7_ = —0.5. The 
equations of the dashed lines are (j) — tt — arctan(l/r) and p = |7_|\/1 + r^/(l + 
7+r) (see the text). Numerical values of the parameters are taken to be 0(0) = 
7r/4, pg = 2, Pp(0) = 0.1 and p{0) = 1. ^^(0) is successively equal to -10, -2.5, 
0, 2.5 and 10 for the different extremals. 



4 Physical conclusions 

We give some qualitative conclusions on the time-optimal control of two-level 
dissipative systems. The discussion concerns the role of dissipation which can be 
beneficial or not for the dynamics and the robustness with respect to dissipative 
parameters of the optimal control. 

The dissipation effect is well summarized by Fig. [5] In this case F > 7+ 
and one sees that as long as the purity of the state decreases (for < z < 1), 
it is advantageous to use a control field, the dissipation being undesirable. On 
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Figure 13: Plot of the optimal control fields vi (solid line) and V2 (dashed line) 
as a function of time t for the extremal trajectory of Fig. dHwith ^^(0) = 5. 
The equation of the horizontal solid line is w = 0. 




t e 

Figure 14: Same as Fig. [T^but for F = 2.5. The equation of the dashed line is 

P = I7-I/7+- 

the contrary, when the purity starts increasing (for — 7_/7-|_ < z < 0) then the 
dissipation alone becomes more efficient and its role positive. The quickest way 
to accelerate the purification of the state consists in letting the dissipation acts. 
This constitutes a non-intuitive physical conclusion which, however, crucially 
depends on the respective values of F and 7+. For instance, if 7-1- > F then all 
the preceding conclusions are modified. 

The robustness of the optimal control with respect to dissipative parameters 
is illustrated by the double-input control. We give different examples. If 7_ = 
then the integrability of the Hamiltonian and the geometrical properties of the 
extremals are preserved when |F — 7_|.| < 2. If 7_ 7^ then the asymptotic 
behavior of the extremals slightly depends on the parameters F, 7_|_ and 7_ (see 
Propositions m and E]). Fig. [THl and [H] show that the extremal control fields 
have also asymptotic behaviors independent of the dissipation. In case (a), 
the control fields tend to a constant whereas a bang-bang structure appears in 
case (b). This conclusion could be interesting for practical applications where 
robustness with respect to physical parameters and simple control fields are 
needed. In addition, due to the simple structure of the time-optimal synthesis, 
shooting techniques will be particularly efficient to determine the control fields 
especially in case (a). 
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Figure 15: (top) Same as Fig. [12] but for the extremal of Fig. [T3]with ^^(0) — 
2.5. (bottom) Evolution of for the same extremal as a function of t. 
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Figure 16: Plot of the extremals of Fig. [T3] up to the first conjugate point. 
The coordinates 6 of the conjugate points are respectively 3.149, 3.116, 3.332, 
3.386 and 3.535 for ^^(O) equal to -10, -2.5, 0, 2.5 and 10. The equations of the 
horizontal and vertical solid lines are respectively (j) =^ tt/2 and 6 = tt. 
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